home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / POLYSPTP.ZIP / TLEP.FOR < prev   
Text File  |  1985-11-29  |  3KB  |  92 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE TLEP
  5. C
  6. C        PURPOSE
  7. C           A SERIES EXPANSION IN LEGENDRE POLYNOMIALS WITH INDEPENDENT
  8. C           VARIABLE X IS TRANSFORMED TO A POLYNOMIAL WITH INDEPENDENT
  9. C           VARIABLE Z, WHERE X=A*Z+B
  10. C
  11. C        USAGE
  12. C           CALL TLEP(A,B,POL,N,C,WORK)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           A     - FACTOR OF LINEAR TERM IN GIVEN LINEAR TRANSFORMATION
  16. C           B     - CONSTANT TERM IN GIVEN LINEAR TRANSFORMATION
  17. C           POL   - COEFFICIENT VECTOR OF POLYNOMIAL (RESULTANT VALUE)
  18. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  19. C           N     - DIMENSION OF COEFFICIENT VECTORS POL AND C
  20. C           C     - GIVEN COEFFICIENT VECTOR OF EXPANSION
  21. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  22. C                   POL AND C MAY BE IDENTICALLY LOCATED
  23. C           WORK  - WORKING STORAGE OF DIMENSION 2*N
  24. C
  25. C        REMARKS
  26. C           COEFFICIENT VECTOR C REMAINS UNCHANGED IF NOT COINCIDING
  27. C           WITH COEFFICIENT VECTOR POL.
  28. C           OPERATION IS BYPASSED IN CASE N LESS THAN 1.
  29. C           THE LINEAR TRANSFORMATION X=A*Z+B OR Z=(1/A)(X-B) TRANSFORMS
  30. C           THE RANGE (-1,+1) IN X TO THE RANGE (ZL,ZR) IN Z, WHERE
  31. C           ZL=-(1+B)/A AND ZR=(1-B)/A.
  32. C           FOR GIVEN ZL, ZR WE HAVE A=2/(ZR-ZL) AND B=-(ZR+ZL)/(ZR-ZL)
  33. C
  34. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  35. C           NONE
  36. C
  37. C        METHOD
  38. C           THE TRANSFORMATION IS BASED ON THE RECURRENCE EQUATION
  39. C           FOR LEGENDRE POLYNOMIALS P(N,X)
  40. C           P(N+1,X)=2*X*P(N,X)-P(N-1,X)-(X*P(N,X)-P(N-1,X))/(N+1),
  41. C           WHERE THE FIRST TERM IN BRACKETS IS THE INDEX,
  42. C           THE SECOND IS THE ARGUMENT.
  43. C           STARTING VALUES ARE P(0,X)=1, P(1,X)=X.
  44. C           THE TRANSFORMATION IS IMPLICITLY DEFINED BY MEANS OF
  45. C           X=A*Z+B TOGETHER WITH
  46. C           SUM(POL(I)*Z**(I-1), SUMMED OVER I FROM 1 TO N)
  47. C           =SUM(C(I)*P(I-1,X), SUMMED OVER I FROM 1 TO N).
  48. C
  49. C     ..................................................................
  50. C
  51.       SUBROUTINE TLEP(A,B,POL,N,C,WORK)
  52. C
  53.       DIMENSION POL(1),C(1),WORK(1)
  54. C
  55. C        TEST OF DIMENSION
  56.       IF(N-1)2,1,3
  57. C
  58. C        DIMENSION LESS THAN 2
  59.     1 POL(1)=C(1)
  60.     2 RETURN
  61. C
  62.     3 POL(1)=C(1)+B*C(2)
  63.       POL(2)=A*C(2)
  64.       IF(N-2)2,2,4
  65. C
  66. C        INITIALIZATION
  67.     4 WORK(1)=1.
  68.       WORK(2)=B
  69.       WORK(3)=0.
  70.       WORK(4)=A
  71.       FI=1.
  72. C
  73. C        CALCULATE COEFFICIENT VECTOR OF NEXT LEGENDRE POLYNOMIAL
  74. C        AND ADD MULTIPLE OF THIS VECTOR TO POLYNOMIAL POL
  75.       DO 6 J=3,N
  76.       FI=FI+1.
  77.       Q=1./FI-1.
  78.       Q1=1.-Q
  79.       P=0.
  80. C
  81.       DO 5 K=2,J
  82.       H=(A*P+B*WORK(2*K-2))*Q1+Q*WORK(2*K-3)
  83.       P=WORK(2*K-2)
  84.       WORK(2*K-2)=H
  85.       WORK(2*K-3)=P
  86.     5 POL(K-1)=POL(K-1)+H*C(J)
  87.       WORK(2*J-1)=0.
  88.       WORK(2*J)=A*P*Q1
  89.     6 POL(J)=C(J)*WORK(2*J)
  90.       RETURN
  91.       END
  92.